library(ggplot2)
library(ggthemes)

power_calculator <- function(mu_t, mu_c, sigma, alpha, N){ 
  lowertail <- (abs(mu_t - mu_c)*sqrt(N))/(2*sigma) 
  uppertail <- -1*lowertail 
  beta <- pnorm(lowertail- qnorm(1-alpha/2), lower.tail=TRUE) + 1- pnorm(uppertail- qnorm(1-alpha/2), lower.tail=FALSE) 
  return(beta) 
} 

#The following were obtained using a pilot study
mu_c <- 61
sigma <- 23
alpha <- 0.05
N <- 750 #Two groups: treated and control, each with size 375 (the rough size of each group given N = 1500, 1 control group, and 3 treatments)

#NOTE: ALL DVS USE THE FEELING THERMOMETER

power <- c()

j <- 1
for(i in seq(0.1, 10, 0.1)){
  power[j] <- power_calculator(mu_t = mu_c + i, mu_c = mu_c, sigma = sigma, alpha = alpha, N = N)
  j <- j + 1
}

df <- data.frame(seq(0.1, 10, 0.1), power, stringsAsFactors = F)
names(df) <- c("effect_size", "power")

cutoff <- data.frame( x = c(-Inf, Inf), y = 0.8, cutoff = factor(50) )

ggplot(data = df) + #geom_line(aes(x = effect_size, y = power1, color = "Including Benefits, N = 1500")) + 
  geom_line(aes(x = effect_size, y = power, color = "a")) + 
  labs(title = "Power Analysis", x = "Effect Size", y = "Power") + theme_minimal() + 
  geom_line(aes( x, y), cutoff) +
  theme_bw() + 
  scale_colour_hc() +
  theme(legend.position = "none")
ggsave("figures/figureC12.pdf", height = 4, width = 4)